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Abstract 

This  paper  explores  the  use  of  a  bio-inspired  control  algorithm,  termed  “local  pursuit”,  as  a  numerical  tool  for  computing 
optimal  control-trajectory  pairs  in  settings  where  analytical  solutions  are  difficult  to  obtain.  Inspired  by  the  foraging  activities 
of  ant  colonies,  local  pursuit  has  been  the  focus  of  recent  work  on  cooperative  optimization.  It  allows  a  group  of  agents  to  solve 
a  broad  class  of  optimal  control  problems  (including  fixed  final  time,  partially-constrained  final  state  problems)  and  affords 
certain  benefits  with  respect  to  the  amount  of  information  (description  of  the  environment,  coordinate  systems,  etc.)  required 
to  solve  the  problem.  Here,  we  present  a  numerical  optimization  method  that  combines  local  pursuit  with  the  well-known 
technique  of  multiple  shooting,  and  compare  the  computational  efficiency  and  capabilities  of  the  two  approaches.  The  proposed 
method  method  can  overcome  some  important  limitations  of  multiple  shooting  by  solving  an  optimal  control  problem  “in 
small  pieces” .  Specifically,  the  use  of  local  pursuit  increases  the  size  of  the  problem  that  can  be  handled  under  a  fixed  set  of 
computational  resources.  Furthermore,  local  pursuit  can  be  effective  in  some  situations  where  multiple  shooting  leads  to  an  ill- 
conditioned  nonlinear  programming  problem.  The  trade-off  is  an  increase  in  computation  time.  We  compare  our  pursuit-based 
method  with  direct  multiple  shooting  using  an  example  that  involves  optimal  orbit  transfer  of  a  simple  satellite. 


Key  words:  Co-operative  control,  Optimization,  Agents,  Group  work,  Trajectories,  Numerical  methods,  Nonlinear 
programming 


1  Introduction 

Over  the  past  decade,  researchers  have  been  increasingly 
looking  to  cooperative  strategies  as  a  means  for  address¬ 
ing  systems  problems  which  are  difficult  to  solve  by  sin¬ 
gle  systems  [18,5] .  The  push  for  cooperation  is  partly  due 
to  the  maturation  of  technology  that  makes  it  possible 
to  develop  meaningful  cooperation  among,  say,  groups  of 
wheeled  or  flying  vehicles,  and  partly  because  of  the  obvi¬ 
ous  success  of  biological  cooperative  systems  which  have 
apparently  evolved  to  “perform  as  more  than  the  sum 
of  their  parts”  [6,13].  Notable  examples  include  worker 


*  This  work  was  supported  by  the  National  Science  Founda¬ 
tion  under  Grant  No.  EIA0088081  and  by  ARO  ODDR&E 
MURI01  Grant  No.  DAAD19-01-1-0465,  (Center  for  Com¬ 
municating  Networked  Control  Systems,  through  Boston 
University). 

*  Corresponding  author.  Tel:  +30-2310-891721,  Fax:  +30- 
2310-891290. 

Email  addresses:  cshao@glue.umd.edu  (Cheng  Shao), 
dcv@uom.gr  (Dimitrios  Hristu-Varsakelis). 


honey  bees,  which  share  information  by  “dancing”  and 
distribute  themselves  among  different  flowers  according 
to  the  “profitability”  of  each  location;  schools  of  fish 
that  swim  in  agile,  tight  formations;  and  ants,  which  use 
pheromone  secretions  for  recruiting  nest-mates  and  for 
discovering  efficient  paths  between  their  nest  and  food 
sources  [4] .  Observations  of  these  and  other  natural  col¬ 
lectives  have  motivated  a  series  of  works  on  the  model¬ 
ing  of  movement  in  animal  groups[4,2,9]  as  well  as  other 
work  on  cooperative  control  strategies,  from  distributed 
collective  covering  and  searching  [18,14],  to  estimating 
by  groups  [15,11],  and  biologically-motivated  optimiza¬ 
tion  [5,3,8]. 

In  particular,  [2]  presented  a  decentralized  organizing 
principle,  inspired  by  the  foraging  behaviors  of  ant 
colonies,  that  allowed  a  group  to  optimize  an  initial 
path  between  two  locations  on  M2 ,  by  having  a  sequence 
of  group  members  travel  from  one  location  to  the  other 
while  each  member  points  its  velocity  vector  towards 
its  predecessor.  This  so-called  “local  pursuit”  rule  was 
generalized  to  non-Euclidean  environments  [7],  and 
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later  to  a  pair  of  pursuit-based  algorithms  that  applied 
to  a  much  broader  class  of  optimal  control  problems 
and  to  systems  with  non-trivial  dynamics [8, 16].  These 
algorithms  require  each  agent  in  the  collective  to  solve 
a  series  of  locally  optimal  control  problems  within  small 
neighborhoods  that  are  defined  by  its  current  state  and 
the  state  of  its  predecessor.  Doing  so  defines  a  con¬ 
trol/trajectory  sequence  that,  under  certain  conditions, 
converges  to  the  optimum. 

Local  pursuit  was  originally  conceived  as  a  means  of 
solving  optimal  control  problems  in  settings  where 
mapping,  communication  and  sensing  capabilities  were 
severely  limited.  The  algorithm  manages  to  avoid  the 
need  for  global  information  by  breaking  up  the  prob¬ 
lem  into  many  pieces,  each  to  be  optimized  by  leader- 
follower  agent  pairs.  Its  key  feature  is  a  reduction  in  the 
range  (measured  by  time,  distance  or  other  metric)  over 
which  computations  must  take  place,  by  paying  a  price 
in  terms  of  the  number  of  agents  that  are  necessary  to 
carry  out  the  algorithm.  For  example,  in  order  for  the 
collective  to  solve  an  optimal  control  problem,  it  is  not 
necessary  to  have  available  a  map  of  the  environment, 
an  agreed-upon  coordinate  system,  or  even  the  coordi¬ 
nates  of  the  target  state;  only  an  initial  feasible  control 
is  required. 

Up  to  now,  discussions  of  local  pursuit  have  focused  on 
broadening  the  domain  of  applicability  of  the  algorithm 
and  on  the  limiting  properties  of  the  agents’  trajecto¬ 
ries.  However,  its  limited  information  requirements  make 
it  a  potentially  useful  tool  in  numerical  trajectory  op¬ 
timization,  where  there  are  often  similar  trade-offs  to 
be  made.  In  particular,  given  a  control  system  which 
must  be  steered  between  two  states  in  a  fixed  amount  of 
time,  one  typically  proceeds  by  “sampling”  the  system 
trajectory  at  pre-determined  times,  and  solving  a  non¬ 
linear  programming  (NLP)  problem  to  determine  the 
best  placement  of  the  samples  in  the  statespace,  subject 
to  constraints  given  by  the  equations  of  motion  and  con¬ 
tinuity,  among  others.  In  that  setting,  there  is  a  balance 
that  must  be  struck  between  the  number  of  trajectory 
samples  and  the  size  of  the  time  intervals  that  separate 
them.  On  one  hand,  we  would  like  to  keep  the  number  of 
segments  small,  so  that  the  associated  NLP  problem  has 
reasonable  storage  requirements.  If  the  number  of  seg¬ 
ments  is  too  small  however,  then  propagating  the  state 
vector  from  one  sample  point  to  the  next  may  require 
high-order  approximation  of  the  equations  of  motion, 
and  thus  become  time-consuming.  If,  on  the  other  hand, 
the  trajectory  is  sampled  densely,  then  there  may  not 
be  adequate  memory  to  solve  the  associated  nonlinear 
programming  (NLP)  problem,  which  at  the  same  time 
is  prone  to  a  more  significant  accumulation  of  numerical 
errors. 

The  purpose  of  this  paper  is  to  explore  the  application 
of  local  pursuit  as  a  computational  tool  in  order  to  ad¬ 
dress  the  shortcomings  outlined  above.  We  introduce 


a  numerical  optimization  method,  titled  “pursuit-based 
multiple  shooting”  (PBMS),  that  combines  a  recently 
reported  local  pursuit  strategy  with  established  numer¬ 
ical  methods.  In  particular,  we  propose  a  modified  ver¬ 
sion  of  the  well-known  multiple-shooting  (MS)  method 
that  uses  local  pursuit  to  overcome  some  of  the  limita¬ 
tions  of  the  original  MS  technique.  The  PBMS  method 
that  is  proposed  here  involves  simulating  group  of  agents 
which  cooperate  to  solve  (numerically)  an  optimal  con¬ 
trol  problem.  Agents  use  standard  MS  to  optimize  their 
short-term  behavior. 

A  comparison  of  the  computational  complexity  and  size 
of  the  resulting  problems  formulated  via  PBMS  ver¬ 
sus  standard  MS,  reveals  that  cooperation  among  group 
members  can  overcome  some  important  limitations  of 
MS.  Our  pursuit-based  formulation  allows  one  to  always 
operate  in  the  domain  of  manageable-sized  problems, 
while  still  being  able  to  treat  problems  with  very  large 
numbers  of  segments.  As  a  result,  it  enlarges  the  maxi¬ 
mum  size  of  problems  that  can  be  handled  given  a  fixed 
amount  of  storage  space.  In  addition,  it  can  reduce  com¬ 
putational  errors  due  to  ill-conditioning  of  the  non-linear 
programming  problem  that  must  be  solved  when  MS  is 
used.  The  trade-off  is  that,  for  well-conditioned  prob¬ 
lems,  demands  more  running  time  than  MS  in  order  to 
converge. 

The  remainder  of  this  paper  is  organized  as  follows:  Sec¬ 
tion  2  defines  a  local  pursuit  algorithm  which  is  applica¬ 
ble  to  optimal  control  problems  with  fixed  final  time  and 
partially-constrained  final  states,  and  gives  the  main  re¬ 
sults  concerning  the  behavior  of  a  collective  under  that 
algorithm.  Section  3  begins  with  a  brief  introduction  to 
numerical  optimization  and  goes  on  to  discuss  the  poten¬ 
tial  advantages  of  modifying  multiple  shooting  (one  of 
the  best-known  and  widely  used  methods  for  numerical 
optimal  control)  to  include  local  pursuit.  An  illustrative 
example  is  presented  in  Section  4  where  we  compare  the 
performance  of  multiple  shooting  to  that  of  puruit-based 
multiple  shooting  in  a  satellite  orbit  transfer  problem. 

2  A  Bio-inspired  Algorithm  for  Optimal  Con¬ 
trol 

Local  pursuit  begins  with  a  group  of  cooperating 
“agents” ,  where  the  term  “agent”  refers  to  a  copy  of  a 
dynamical  system: 

Xk  =  f(xk,uk ),  xk(t)  g  Rn,uk(t)  g  n  C  Mm  (1) 

for  k  =  0, 1, 2  . . ..  The  problem  of  interest  is  as  follows: 

Problem  1  Find  a  trajectory  x*(t),  and  a  final  state 
x*  (T) ,  (T  fixed)  that  minimize 

rto+T 

J(x,x,t0,T)  =  /  g(x,x)dt  + G(x(t0  +  T))  (2) 

d  to 


subject  to  the  constraints  a; (to)  =  %o,  and  Q(x(to+T ))  = 

0. 

Here  it  is  assumed  that  g (x (t),  x (t))  >  0,  G(x(to+T ))  > 
0  and  that  Q(-)  is  an  algebraic  function  of  the  state. 

Definition  1  Given  the  final  state  constraint  Q(x)  =  0, 
the  constraint  set  of  x  is 

Sq  =  {x  &  M"|<2(x)  =  0}. 

For  any  pair  of  fixed  states  a,  b  £  D  C  R",  suppose 
the  optimal  trajectory  from  a  to  b  with  fixed  final  time 
(minimizing  J  with  respect  to  x  only)  is  denoted  by 
x*(t).  Then  the  cost  of  following  x*(-)  is  denoted  by: 

v(a,b,t0,T)=  [  g(x*(t),x*(t))dt  (3) 

Jt0 

subject  to  x(to)  =  a ,  x(to  +  T)  =  b. 

Now,  let  x*(t)  be  the  optimal  trajectory  (over  T  units 
of  time)  from  an  initial  state  a  to  the  constraint  set  Sq. 
The  cost  of  following  x*(-)  is  denoted  by: 

nto~\-T 

rjQ(a,t0,T)=  /  g(x*,x*)dt  + G{x*{t0  +  T)) 

Jt0 

=  minJ(x,x,to,T)  (4) 

X 


subject  to  x(t0)  =  a,  Q(x(t0  +  T))  =  0. 

2.1  Algorithm 

In  previous  works  [8,16]  we  have  proposed  two  classes 
of  algorithms,  one  for  problems  with  fixed  boundary 
conditions  (final  time  and  final  states),  and  another  for 
problems  with  partially-constrained  boundary  condi¬ 
tions.  We  will  briefly  present  a  “sampled”  version  local 
pursuit  for  trajectory  optimization  problems  with  fixed 
final  time  and  partially-constrained  final  states  (based 
on  algorithms  from  [8,16]). 

We  assume  that  there  is  an  available  initial  (but  sub- 
optimal)  feasible  control/trajectory  pair  (ufeas(t),Xfeas(t)) 
for  (1),  which  could  have  been  obtained  through  explo¬ 
ration  or  from  a-priori  knowledge.  Agents  are  scheduled 
to  leave  the  starting  state  x0  sequentially,  separated  by 
a  pursuit  interval  of  A  time  units,  i.e. ,  if  the  first  (ini¬ 
tial)  agent  leaves  at  time  t o,  the  kth  agent  will  leave  at 
tk  =  t o  +  kA,  k  =  0,1,2,....  The  first  agent  is  required 
to  follow  x feas  from  xq  to  Sq  .  The  next  agent  attempts 
to  intercept  its  predecessor  along  an  optimal  trajecto¬ 
ries  defined  by  (3),  as  long  as  the  predecessor  has  not 
yet  reached  the  target  set  Sq.  If  the  predecessor  is  al¬ 
ready  in  Sq  then  the  follower  moves  along  the  optimal 


trajectory  defined  by  (4)  (see  Fig.  1  for  illustration). 
Agents  do  not  monitor  their  predecessors  continuously, 
but  instead  updates  their  trajectories  with  the  sampling 
rate  of  1/6.  The  precise  rules  that  govern  the  movement 
of  each  agent  are: 

Algorithm  1  (Sampled  Local  Pursuit):  Identify  the 
starting  state  x$  on  D  and  the  constraint  set  Sq. 
Let  xo (t)  ( t  £  [0,To])  be  an  initial  trajectory  satisfy¬ 
ing  (1)  with  cco(0)  =  Xo,  <3(a;o(Tb))  =  0.  Choose  the 
pursuit  interval  A  and  updating  interval  S  such  that 
0  <  8  <  A  <  T0. 

(1)  For  k  =  1,  2,  3  . . let  tk  =  fcA  be  the  starting  time 
of  the  kth  agent,  i.e.  uk(t)  =  0,  xk{t)  =  xq  for 
0  <t<tk. 

(2)  When  t  =  tk  +  iS,i  =  0,1, 2,  3,...,  calculate  the 
control  u)(t)  that  achieves  (subj.  to  (1)): 

v(xk(t),Xk-i(t),t,  A),  if  xk-i(t)  i  SQ 

<  (t  €  [M  + A]) 

V Q(Xk(t),t,T  -  ( t  -  tk)),  if  Xk-i(t)  e  Sq 

(t  £  [ t ,  tk  +  T]) 

(3)  Apply  Uk(t)  =  u(k+i6(t)  to  the  kth  agent  for  t.  £ 
[tk  +  iS,tk  +  {i  +  1)5)  if  A  +  iS  <  T,  or  for  t  £ 
[tk  +  id,  tk  +  T)  otherwise. 

(4)  Repeat  from  step  2  until  the  kth  agent  reaches  Sq. 


Initial  Trajectory 


Fig.  1.  Illustration  of  local  pursuit  on  a  manifold.  Each  agent 
except  the  first  one  is  trying  to  “catch  up”  its  predecessor 
before  the  predecessor  reaches  Sq,  and  trying  to  reach  Sq 
via  a  locally  optimal  trajectory  otherwise. 

The  two  adjustable  parameters  in  sampled  local  pursuit 
(SLP )  algorithm  are  the  pursuit  interval  A ,  which  defines 
how  frequent  agents  depart  from  xq  (consequently  how 
far  away  agents  are  separated) ,  and  the  updating  interval 
S,  which  defines  the  frequency  with  which  each  agent 
updates  its  trajectory.  We  will  refer  to  two  successive 
agents  as  a  “pursuit  pair”.  Within  a  pair,  the  (k  —  l)th 
agent  will  be  known  as  the  “leader”,  and  the  kth  agent 
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as  the  “follower”.  We  will  refer  to  the  times  t\  =  tk  + 
id,  i  =  0, 1,2,3...  as  the  “updating  times”.  The  SLP 
algorithm  employs  two  types  of  pursuit  -  “catching  up” 
and  “free  running”  -  depending  on  whether  the  leader 
has  reached  the  final  constraint  set  Sq  or  not.  The  former 
lets  agents  “learn”  from  their  leaders,  meanwhile  the 
latter  enables  them  to  find  the  optimal  final  state.  This  is 
different  from  the  version  of  sampled  local  pursuit  in  [8], 
which  involves  only  “catching  up”  stages,  and  because 
5  is  fixed,  rather  than  determined  on-line  (as  is  the  case 
in  [16]),  the  total  number  of  updates  performed  by  each 
agent  is  pre-determined. 

We  now  state  the  main  result  concerning  the  limiting 
behavior  of  the  SLP  algorithm.  Before  proceeding  to 
the  main  theorem,  we  must  require  that  the  cost  of  (2) 
changes  “little”  for  small  changes  to  the  endpoints  of  a 
trajectory: 

Condition  1  Assume  for  a  generic  trajectory  X\ (i) 
there  exists  an  e  >  0  such  that  for  all  a,  b±,  62  G  D  and  all 
A  >  0,  there  exists  a  trajectory  X2 (t)  such  that  the  cost 
C(xi,0,T)  (2:1(0)  =  a,x±(T)  =  61)  from  a  to  b\  and  cost 
C(x2,0,T)  (2:2(0)  =  a,X2(T)  =  62)  from  a  to  62  satisfy 

\\bi  -  &2II00  <  £ 

=*  \\C(Xl,0,T)-  C(x2, 0,  T)  IU  <CA 

for  some  constant  C  independent  of  A. 


3.1  A  Brief  Review  of  Multiple  Shooting 

In  principle,  MS  is  a  type  of  NLP  method [1].  Suppose, 
for  now,  that  we  are  interested  in  minimizing  the  scalar 
cost  function  F(x)  :  — >  K,  subject  to  m  (m  <  n) 

equality  constraints  c(x)  =  0,  where  c(x )  is  an  (m  x  1) 
vector.  To  do  this,  we  first  introduce  the  Lagrangian: 

L(x,  A)  =  F(x)  +  Atc(x)  (5) 

The  necessary  conditions  for  the  point  (x*,  A*)  to  be  an 
optimum  are  satisfied  by  the  stationary  points  of  the 
Lagrangian  1 : 

\/xL(x,  A)  =  A7xF(x)  +  Vxc(x)TX  =  0  (6) 

and 

V\L(x,  A)  =  c(x)  =  0  (7) 

The  equations  (6)~(7)  are  usually  solved  via  Newton’s 
method.  Proceeding  formally  we  obtain  the  following 
Karush-Kuhn-Tucker  system: 


Hl 

Vxc(x)T 

Ax 

-VxF(x)  -  Vxc( x)TA 

_Vxc(x) 

0 

_  AA_ 

-c{x) 

(8) 


Theorem  1  Suppose  a  group  of  agents  evolve  under 
sampled  local  pursuit,  and  that  at  every  updating  time 
tlk,  the  locally  optimal  trajectories  from  follower  to  leader 
are  unique.  If  the  updating  interval  S  and  pursuit  inter¬ 
val  A  satisfy  0  <  6  <  A  and  Condition  1  holds,  then  the 
limiting  trajectory  obtained  is  unique  and  locally  opti¬ 
mal.  It  is  also  smooth  if  the  locally  optimal  trajectories 
calculated  at  every  updating  time  are  smooth. 


where  Ax  is  the  “search-direction”  and  Hl  is  the  Hessian 
of  the  Lagrangian 

m 

HL  =  VlF(  x)  +  J2^lci  (9) 

i=  1 

For  convenience,  we  define 


PROOF.  See  [8,16], 


Hl  Vxc(  x)T 
\7xc(x)  0 


(10) 


3  Local  Pursuit  as  a  Numerical  Computational 
Tool 

Multiple-shooting,  together  with  its  various  improved 
forms,  could  be  considered  a  workhorse  of  numerical  op¬ 
timization.  However,  the  large  storage  requirements  as¬ 
sociated  with  MS  limit  the  so-called  “permitting  size”  of 
problems  which  can  be  solved  on  a  digital  computer  [1]. 
In  the  following,  we  briefly  review  the  technique  known 
as  multiple  shooting  (MS)  before  showing  how  it  im¬ 
proves  when  combined  with  SLP. 


Suppose  now  that  we  are  interested  in  optimizing  (2) 
subject  to  the  dynamics 


i  =  f(x(t),u{t))- 

1  Here,  we  use  the  notation: 


A7XA(  x) 


da  1 

da  1 

da  1 

dxi 

dx  2 

dxn 

da2 

da2 

da2 

dxi 

dx  2 

dx  n 

dan 

dan 

dan 

dx  1 

dx  2 

dxn 

(11) 
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Doing  so  with  MS  (here  we  mainly  refer  to  di¬ 
rect  MS)  involves  breaking  up  the  trajectory  into 
“shorter”  pieces  [1]  by  partitioning  the  time  domain 
[to,tf]  =  U?=r11[tl,ti+1),  to  =  h  <  ■■  ■  <  tN  =  tf.  Each 
subinterval  [f,;,  tl+-[)  is  called  a  “segment”. 

Multiple  shooting  uses  NLP  to  find  the  optimal  trajec¬ 
tory  (i.e.  to  minimize  a  scalar  function  along  a  trajec¬ 
tory),  by  defining  the  NLP  variables  v  =  [zq, . . . ,  z/jv] 
to  be  the  arguments  -  usually  they  are  concatenations 
of  the  states  and  the  corresponding  controls  -  at  times 
along  a  trajectory.  The  stationary  points 
, . . . ,  obtained  via  NLP  will  approach  points  on 
the  optimal  trajectory.  Without  loss  of  generality,  we 
choose  to  sample  the  state  and  control  vectors  at  the 
segment  endpoints,  and  define  the  NLP  variables 

v  =  {x1,u1,x2,u2,  . . .  ,xN,uN}  (12) 


the  precise  state  evolution  by  integration.  Second,  using 
linear  approximation  reduces  computational  burden;  on 
the  other  hand,  higher-order  approximations  are  likely 
to  produce  more  precise  results,  but  they  will  also  give 
rise  to  more  complicated  equations  in  (6)~(13)  and  the 
resulting  NLP  problem  will  require  more  time  to  solve. 

One  shortcoming  of  linear  approximation  is  that  the  time 
separation  between  neighboring  points  must  be  limited 
to  a  small  time  step  h.  If  h  is  “too  large”,  then  2j+i 
will  not  be  a  good  estimate  of  the  true  state  of  (11)  as 
it  evolves  from  Xi  to  Xi+ 1  under  Uj ,  and  NLP  will  pro¬ 
duce  erroneous  results.  The  estimation  error  is  of  o(hn), 
where  n  £  N  is  the  order  of  the  approximation  methods 
[10].  It  is  clear  that  a  smaller  h  results  in  better  approxi¬ 
mation.  For  that  reason,  keeping  the  segment  size  small 
is  desirable  and  is  a  key  feature  of  MS  compared  with 
other  single-step  shooting  methods  [1] . 


Notice  that  the  dimension  of  Xi  and  Ui,  i  =  1, N  is 
Mx  and  Mu,  respectively. 

As  the  NLP  variables  are  adjusted,  one  must  ensure  that 
they  satisfy  the  the  system  dynamics  (11),  and  that  se¬ 
quential  trajectory  segments  obey  a  matching  condition 
at  their  boundaries.  This  means  that  the  evolution  of 
(11)  from  Xi  at  U  with  input  Ui ,  should  steer  the  system 
to  Xi+ 1  at  ti+ 1.  This  suggests  that  the  following  con¬ 
straints  are  necessary: 


£2  —  2i 

23  -  22 


C(2) 


2 N  -  2jv— 1 


=  0 


</>o(a:l!  to) 
4>t{xn ,  tf) 


(13) 


where  the  functions  4>o(xi,  to)  and  4>T{xN,tf)  represent 
the  initial  and  final  condition  constraints,  respectively, 
and  2 i  are  to  be  calculated  by  integrating  the  differential 
equation  (11)  from  ti  to  ti+i. 

In  practice,  one  often  computes  x i  only  approximately, 
instead  of  integrating  the  complete  equations  of  motion. 
For  example,  Euler’s  method  provides  a  first-order  ap¬ 
proximation  to  the  integration  of  (11): 


xi+i  =2  i  +  hf(Xi,Ui) 


(14) 


Let  us  assume  that  we  can  fix  an  acceptable  h ,  i.e., 
one  that  is  small  enough  to  lead  to  convergence  for  the 
associated  NLP.  Then  the  permitting  size  of  problems 
that  can  be  solved  by  MS  depends  solely  on  the  number 
of  steps  N  required  to  cover  the  time  span  [0,  T],  with 
T  =  (N  —  1  )h.  Thus,  solving  large-scale  problems  (with 
large  time  spans),  requires  a  large  number  of  segments 
N  =  Th  +  1.  For  problems  with  varying  dynamics  but 
fixed  time  span,  the  requirement  for  N  could  be  large 
because  nonlinearities  in  the  dynamics  make  it  neces¬ 
sary  to  use  small  time  steps  to  ensure  that  the  estimates 
x i  are  precise  enough  and  that  the  associated  NLP  will 
convergence. 

Remark  2  For  any  given  set  of  dynamics  and  choice 
of  NLP  method,  the  time  step  h  is  upper  bounded.  The 
permitting  size  in  MS  is  proportional  to  the  number  of 
segments  N,  if  h  is  pre- determined. 

Suppose  now  that  the  dimensions  of  system  state  x  and 
control  u  are  Mx  and  Mu  respectively,  and  the  number 
of  segments  is  N,  then  the  dimension  of  NLP  variables 
for  a  multiple  shooting  method  is  nx  =  ( Mx  +  MU)N, 
and  the  NLP  has  at  least  MXN  constraints  (the  number 
of  constraints  could  be  larger  if  the  final  states  are  fixed) . 
The  Hessian  in  (9)  is  of  dimension  nx  x  nx  =  (Mx  + 
MU)2N2.  To  proceed  with  NLP,  one  must  obtain  the 
solution  of  (8),  which  requires  finding  the  solution  to 
a  linear  system  of  algebraic  equations  with  dimension 
at  least  (Mu  +  2MX)2N2.  Here  we  have  seen  that  the 
number  of  segments  involved  in  the  calculation  affects 
the  “degree  of  labor-consumption”  with  0(N2). 


where  h  is  the  time  step  of  each  segment.  The  choice 
of  approximation  should  be  based  on  the  following  con¬ 
siderations:  First,  we  have  sampled  the  controls  at  the 
segment  boundaries  (the  samples  are  to  be  optimized) 
instead  of  considering  the  continuous-time  control  on 
the  intervals  therefore,  we  can  not  compute 


3.2  Numerical  Optimization  by  Local  Pursuit 

As  we  have  seen,  the  number  of  segments  used  in  MS 
affects  the  dimension  of  the  associated  NLP  variables 
and  the  computational  complexity  of  the  NLP  problem 
to  be  solved.  For  that  reason,  it  would  be  desirable  to 
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use  fewer  segments  with  the  MS  method.  However,  be¬ 
cause  the  time  step  h  is  upper  bounded  for  a  fixed  set 
of  dynamics  and  NLP  algorithm,  decreasing  the  number 
of  segments  means  that  NLP  process  can  only  deal  with 
trajectories  spanning  shorter  time  intervals.  This  limi¬ 
tation  can  be  circumvented  by  combining  local  pursuit 
with  direct  MS,  to  obtain  what  we  refer  to  as  PBMS. 
Specifically,  one  can  introduce  a  sequence  of  simulated 
agents  that  pursue  each  other  using  the  algorithm  given 
in  Sec.  2.  Each  agent  will  use  MS  to  compute  the  opti¬ 
mal  trajectory  from  its  own  state  to  that  of  its  leader, 
giving  rise  to  a  series  of  smaller  NLP  problems,  whose 
time  span  is  limited  by  the  pursuit  interval  A.  Although 
there  will  be  more  of  these  NLP  problems  to  solve,  their 
lower  dimension  will  make  it  possible  to  handle  larger 
optimization  problems  overall. 


3.2.1  Decreasing  the  size  of  the  NLP  problems  during 
computation 


Fig.  2.  Local  pursuit  could  decrease  the  problem  size  in¬ 
volved  in  calculation  at  every  updating  step.  The  number  of 
variables  calculated  in  multiple  shooting  is  N  =  13,  and  the 
number  of  variables  calculated  by  each  agent  in  local  pursuit 
is  JVa  =  5. 

For  simplicity,  Rx  h  =  T/{N  —  1)  and  select  the  pur¬ 
suit  interval  A  =  {Na  —  1  )h,  Na  £  N,  where  Na  is 
the  number  of  segments  within  A.  Usually  we  will  have 
Na  «  N.  For  convenience,  we  also  choose  the  updat¬ 
ing  interval  in  the  SLP  algorithm  to  be  an  integer  multi¬ 
ple  of  segment  size,  <5  =  Nsh,  Ns  £  N  so  that  the  agents 
are  always  updating  their  trajectories  at  the  times  tj 
which  we  chose  to  sample  the  trajectory  of  (11).  At  each 
updating  step  t  =  id,  each  agent  is  solving  a  NLP  over 
Na  segments,  with  a  time  span  of  {Na  —  l)h  instead  of 
(N  —  1  )h,  as  illustrated  in  Fig.  2.  Because  the  computa¬ 
tional  complexity  of  MS  is  related  to  the  square  of  the 
number  of  segments,  using  Na  <<  N  will  significantly 
decrease  the  computational  burden  for  each  agent.  Ta¬ 
ble  1  shows  the  dimensions  of  the  NLP  vector,  v,  the 
constraints  c{x),  and  the  matrix  H  in  Eq.  (10),  respec¬ 
tively,  for  MS  and  PBMS.  The  dimension  of  the  associ¬ 
ated  Karush-Kuhn- Tucker  system  is  on  the  order  of  N2 


Table  1 

Comparing  the  dimensions  of  the  NLP  problem  variables 
when  using  Multiple  Shooting  (MS)  vs.  Pursuit-based  Mul¬ 
tiple  Shooting  (PBMS). 


MS 

PBMS 

dim(iz) 

(Mx  +  MU)N 

(. Mx  +  Mu)Na 

dim(c(x)) 

MXN 

MxNa 

dim)#) 

{{Mu  +  2MX)N)2 

{{Mu  +  2Mx)Na)2 

for  MS,  vs.  N \  for  PBMS.  Operating  on  large  matrices 
consumes  large  amounts  of  memory,  especially  when  us¬ 
ing  non-iterative  methods,  e.g.,  the  memory  of  a  typical 
desktop  PC  can  be  easily  used  up  by  a  matrix  with  di¬ 
mension  of  5000  x  5000  in  Matlab  when  using  64  bits 
of  digital  precision.  However,  under  PBMS,  the  matrix 
size  is  scaled  down  by  a  factor  of  ( Na/N )  x  (Na/N)  and 
hardware  requirements  can  be  decreased  significantly  for 
any  given  problem. 

We  note  that  under  PBMS,  each  agent  needs  to  solve 
(. N  —  Na)/Ns  “smaller”  MS  problems  in  order  to  reach 
the  target  set  Sq  from  the  initial  state  xg,  but  the  time 
needed  to  do  so  denoted  by  Ta  -  will  generally  be  less 
than  the  total  iterative  time  in  multiple  shooting,  which 
we  will  denote  by  Tms-  Of  course,  the  total  running 
time  for  PBMS  -  denoted  by  Tpp  -  may  be  greater  than 
Tms  because  local  pursuit  relies  on  multiple  agents  to 
converge  to  the  optimum.  Our  experience  in  prior  work 
on  local  pursuit  and  in  various  numerical  experiments 
(including  the  example  in  next  Section)  has  been  that, 
for  well-conditioned  problems,  the  convergence  rate  of 
PBMS  is  usually  slower  than  that  of  MS.  As  expected, 
decreasing  Na  leads  to  slower  convergence  for  PBMS. 
However,  the  added  running  time  comes  with  the  ben¬ 
efit  of  lower  memory  requirements,  allowing  us  to  han¬ 
dle  problems  with  larger  state  vectors  and  longer  time 
horizons.  At  the  same  time,  if  N a  is  decreased  and  the 
available  storage  is  fixed,  one  can  afford  to  also  decrease 
the  segment  size  h.  Doing  so  has  the  effect  of  improving 
the  state  estimates  Xj  used  to  formulate  the  NLP,  with¬ 
out  making  the  problem  ill-conditioned.  This  situation 
will  be  illustrated  in  the  next  Section. 

3.2.2  Reducing  numerical  error  when  H  is  ill- 
conditioned 

Besides  maximum  problem  size,  another  important  con¬ 
sideration  in  numerical  optimal  control,  is  the  computa¬ 
tional  error  introduced  by  the  finite  precision  of  digital 
computers  and  by  algorithmic  accuracy.  It  is  possible 
that  the  error  is  too  large  to  obtain  useful  results,  e.g., 
solving  a  linear  system  with  large  condition  number  can 
result  in  unacceptable  errors  and  non-convergence.  In 
such  settings,  correction  algorithms,  such  as  Tikhonov 
regularization,  can  be  applied  [10,12];  they  are,  however, 
time  and  storage  consuming,  and  do  not  always  succeed. 

For  the  optimal  control  problem  of  interest,  consider,  for 
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example,  using  Gaussian  elimination  to  solve  (8)  with 
limited  digital  precision.  Every  step  of  the  elimination 
algorithm  introduces  some  truncation  error,  so  that  the 
total  accumulated  error  when  solving  a  large  linear  sys¬ 
tem  will  generally  be  much  larger  than  that  associated 
with  solving  one  of  lower  dimension,  because  the  number 
of  steps  required  by  the  algorithm  is  proportional  to  the 
system  size.  Furthermore,  if  H  in  (10)  is  ill-conditioned, 
then  the  numerical  solution  of  (8)  introduces  small  er¬ 
rors  which  accumulate  towards  the  final  segment  N.  If 
the  problem’s  time  horizon  is  long,  the  accumulated  er¬ 
ror  may  lead  to  erroneous  results,  or  prevent  MS  from 
converging.  PBMS  can  help  reduce  these  numerical  er¬ 
rors  because  the  algorithm’s  simulated  agents  solve  MS 
problems  with  a  shorter  time  horizon,  compared  to  that 
of  the  original  problem.  This  implies  that  the  dimension 
of  (8)  for  each  agent  is  reduced  and  there  will  be  cases  in 
which  PBMS  will  succeed  where  MS  failed  to  converge. 
Furthermore,  if  every  locally  optimal  trajectory  satisfies 
the  convergence  criteria  of  the  numerical  method  used 
to  solve  the  “short-range”  MS  problems  between  leader- 
follower  pairs,  then  the  convergence  of  the  agents’  tra¬ 
jectory  sequence  is  guaranteed  by  the  local  pursuit  algo¬ 
rithm  itself. 


4  Example:  An  Orbit  Transfer  Problem 

Consider  an  idealized  spacecraft  which  must  be  trans- 
fered  from  one  stable  orbit 2  to  another,  within  some 
fixed  time  T.  For  simplicity,  we  only  consider  the  effect 
of  the  Earth  and  restrict  the  problem  to  a  plane,  as  il¬ 
lustrated  in  in  Fig.  3. 


Fig.  3.  A  planar  spacecraft  in  orbit  around  the  Earth. 


3.2.3  Remarks 


In  summary,  the  combination  of  MS  with  local  pursuit 
can  increase  the  permitting  size  of  problems  that  can 
be  handled  with  fixed  storage,  because  at  every  updat¬ 
ing  step,  each  agent  deals  with  a  problem  with  “reduced 
size”.  The  development  of  NLP  algorithms  (and  of  MS) 
has  followed  the  growth  of  the  digital  computer.  The  size 
of  a  typical  application  in  the  early  1960s  was  n,  m  sa  10, 
while  in  the  1970s  and  early  1980s  most  application  were 
of  size  as  n,  m  <  100.  With  subsequent  advances  in  lin¬ 
ear  algebra  techniques,  such  as  matrix  sparsity,  and  on¬ 
going  progress  in  the  semiconductor  industry,  the  per¬ 
mitting  size  in  late  1990s  was  n,m  ss  10,000  [1],  How¬ 
ever,  when  using  a  fixed  time  partition,  PBMS  involves 
solving  problems  of  size  N&  instead  of  N.  Therefore,  we 
can  address  much  larger  problems  under  the  limits  im¬ 
posed  by  the  hardware.  Although  it  requires  more  run¬ 
ning  time,  PBMS  does  provide  a  feasible  solution  when 
the  traditional  formulation  exceeds  those  limits,  making 
it  impossible  to  proceed. 


The  dynamics  of  this  system  are 


Pui 

mg 


2  dr  Puo 


mgr 


P{u\+ul) 


U\  =  I  sin(yj) 
u2  =  I  cos(<p) 


(15) 


where  r  is  the  distance  between  the  spacecraft  and  the 
center  of  Earth,  6  is  its  longitude  with  respect  to  the 
horizontal  line,  m  is  the  mass  of  the  spacecraft,  ue  is 
the  gravitational  parameter  of  Earth,  P  is  a  constant 
concerning  engine  power,  and  g  is  the  acceleration  of 
gravity  at  sea  level  [17].  The  control  inputs,  u\  and  u2, 
are  functions  of  /  and  ip,  the  thrust  force  and  thrust 
steering  angle  with  respect  to  the  tangent  of  the  local 
orbit. 


We  would  like  to  minimize  the  fuel  consumed  (equiva¬ 
lently,  to  maximize  the  final  mass  m(T)): 


In  cases  where  MS  fails  to  converge  because  the  er¬ 
rors  introduced  by  the  approximation  to  (11)  are  too 
great,  PBMS  may  succeed  by  reducing  the  segments  N&, 
thereby  reducing  the  accumulated  error  over  the  tra¬ 
jectory  of  a  single  agent.  In  next  section  we  present  an 
example  of  MS  stagnancy  caused  by  an  ill-conditioned 
matrix  H  and  how  PBMS  can  avoid  the  problem. 


J=  [  (ui(f)2  +  u2{t)2)dt 
Jo 


2  The  term  “stable  orbit”  means  that  the  spacecraft  will 

remain  on  this  orbit  if  no  external  force  (other  than  gravity) 
is  applied,  i.e.  9  =  ue/t 3. 
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while  steering  the  system  (15)  from  the  initial  condition 
r(0)  =  Rq,  r(0)  =  0,  0(0)  =  \Jue/Rq,  m(0)  =  M0  to  the 
final  condition  r(T)  =  Rx,f(T)  =  0 ,  0(T)  =  s/ueJR^., 
where  T  is  fixed.  For  simplicity,  we  assumed  there  was 
no  upper  bound  of  the  thrust  force,  thus  there  is  no 
restriction  for  the  control  U\  and  u-i. 

4-1  Applying  MS  vs.  PBMS 


IE  —  8  (to  guarantee  little  improvement  with  future  it¬ 
erations)  and  ||c(x)||  <  (IE  —  15)  x  m,  where  to  was  the 
dimension  of  constraints  (to  guarantee  that  the  trajec¬ 
tory  satisfies  the  system  dynamics).  The  total  time  T 
was  fixed  to  300  minutes.  The  numerical  properties  of 
the  problem  -  and  consequently  the  performance  of  the 
two  methods  -  dependeded  on  the  selection  of  the  total 
number  of  segments  TV. 


To  solve  the  problem  by  MS,  we  used  Euler’s  method 
to  estimate  the  piecewise  integration  of  Eq.  (11)  and 
imposed  the  following  constraints: 


0  =  r3+ 1  -  fj 
=  rj+ 1  -  (jj  +  h3rj) 

0  =  f3+ 1  -  fj 

_  .  (  .  .  ,  M  UE 

—  ri+ 1  (rj  +  /ij(  2 

r3  r3 

o  =  e3+1-~e3 


Puij 

m3g 


)) 


—  0j+i  —  (0j  +  h3(— 


2d3f3 


Pu?j 

m3gr3 


)) 


0  =  rrij+ 1  —  m3 
=  rrij+i  —  ( m3  —  h3 
for  j  =  1,2,3..., 


P(u\  +  ul) 

a2 

N  —  1  and 


) 


o  =  n  -  Ro 

o  =  f  -i  -  ii,{) 

o  =  0i  -  0O 
0  =  mi  —  M0 


0  =  r  ^  —  Rt 
0  =  r  ^  —  Rt 
0  =  6t  —  0t 


(16) 


where  r;  =  r((i  —  1  )/i),fj  =  f((i  —  1  )h), . . .  and  h  = 
T/(N  —  1)  was  the  time  step.  The  dimension  of  con¬ 
straints  c(x )  (to  as  we  denoted  before)  was  (4  x  ( N  — 
1)  +7).  The  NLP  variable 


4-2.1  Well- conditioned  case 

With  N  =  101,  the  matrix  in  Eq.  (8)  was  well  condi¬ 
tioned  (can  be  verified  by  its  condition  number  and  the 
fast  convergence  rate  for  MS  shown  below).  The  per¬ 
formance  of  MS  and  PBMS  are  summarized  in  Table 
2,  where  c (H)  was  the  condition  number  of  matrix  H\ 
The  “iterations”  column  lists  the  iteration  numbers  for 
convergence  in  multiple  shooting,  and  number  of  agents 
needed  for  local  pursuit  to  converge,  respectively. 


Table  2 

Comparison  between  Multiple  Shooting  and  Local  Pursuit 
in  well-conditioned  case 


N=101 

MS 

PBMS 

PBMS 

? 

II 

CO 

o 

(Na  =  60, 

Ns  =  16) 

CA 

CO 

II 

£ 

Comp.  Time 

66.8594 

793.8594 

430.2344 

Iterations 

14 

277 

41 

m(T ) 

0.524246124 

0.524245714 

0.524246100 

l|c(*)|| 

3.7144E-14 

9.2499E-15 

9.9170E-15 

Ave(c(H)) 

1.5931E+6 

4.8919E+5 

5.3612E+5 

We  can  see  that  both  methods  were  successful  and  that 
the  convergence  rate  of  PBMS  was  slower  than  that  of 
MS.  Increasing  Na  resulted  in  increased  need  for  storage 
and  decreased  running  time.  If  Na  =  N,  then  PBMS 
reduces  to  MS.  The  final  trajectories  obtained  by  both 
methods  are  shown  in  Fig.  4. 


V  =  {n , . . . ,  rjv,  h , . . . ,  fN ,  0t ,  ■  •  • ,  0JV ,  TO! , . . . ,  mN 

!  U\i,  .  .  .  ,  Uljv-l,  W21,  ■  ■  •  ,  M2JV-l}  (17) 

was  of  dimension  (6  x  TV  —  2) . 

When  applying  local  pursuit,  the  TV  should  be  replaced 
by  Na  in  the  above  equations.  The  dimension  of  the 
constraint  set  and  the  NLP  variable  were  (4  x  (Na  — 
1)  +  7)  and  (6  x  Na  —  2),  respectively. 

4-2  Results 

We  solved  the  problem  both  by  PBMS  and  standard  MS. 
The  criteria  for  convergence  were  ||m(T)j+i  —  to(T)j||  < 


4-2.2  Ill-conditioned  case 

When  the  segment  size  was  halved,  i.e. ,  TV  >  201,  MS 
became  stagnant  because  the  matrix  H  in  (10)  was  ill- 
conditioned.  The  error  generated  by  solving  Karush- 
Kuhn- Tucker  system  became  so  large  that  the  NLP  al¬ 
gorithm  (using  Newton’s  method)  was  not  able  to  con¬ 
verge.  The  large  condition  number  of  H  (our  criteria  of 
ill-condition)  was  due  to  the  large  number  of  segments 
(in  fact  H's  condition  number  increased  with  the  size  of 
segments,  see  Table  3).  For  PBMS,  the  reduction  in  the 
number  of  segments  for  the  sub-problems  solved  by  pur¬ 
suing  agents  meant  a  reduction  in  the  condition  number 
of  H  when  using  small  Na  ■  The  results  from  both  meth¬ 
ods  are  summarized  in  Table  3. 


Fig.  4.  Trajectories  of  both  methods  under  well-conditioned 
case  with  N  =  101,  Na  =  30,  Ns  =  16.  The  trajectories 
obtained  from  both  methods  were  virtually  identical. 


Table  3 

Comparison  between  Multiple  Shooting  and  Local  Pursuit 
in  ill-conditioned  case _ 


N=201 

MS 

PBMS 

PBMS 

o 

CO 

II 

if. 

(IVa  =  60, 

Ns  =  16) 

cn" 

co 

II 

5? 

Time 

>100000 

32911.8750 

10676.9844 

Iteration 

>3000 

7223 

603 

m(T) 

0.523948746 

0.524560874 

0.524571236 

IK*)II 

0.01709873 

6.0790E-14 

1.4019E-14 

Ave(c  (H)) 

1.8263E+10 

5.0576E+5 

4.7881E+5 

Max(c  (H)) 

4.3377E+12 

6.1142E+5 

5.6897E+5 

N=301 

MS 

PBMS 

PBMS 

o 

CO 

II 

(1VA  =  60, 

to 

i-H 

II 

£ 

Ns  =  32) 

Time 

>360000 

239601.4219 

81264.3906 

Iteration 

>3200 

30722 

2845 

m(T) 

0.527347984 

0.524643889 

0.524700485 

l|c(*)ll 

0.04135162 

2.6051E-13 

1.2922E-13 

Ave(c  (H)) 

4.5106E+10 

5.3032E+5 

4.8184E+5 

Max(c  (H)) 

7.3201E+10 

6.5692E+5 

6.5833E+5 

In  this  case,  MS  could  not  converge  (the  values  of  the 
constraint  residues  c{x)  did  not  become  sufficiently 
small),  and  the  iterative  process  became  stagnant.  On 
the  other  hand,  PBMS  was  effective  in  producing  the 
optimal  trajectory.  The  final  trajectories  of  both  meth¬ 
ods  are  shown  in  Fig.  5.  By  comparing  to  the  trajecto¬ 
ries  generated  in  the  well-conditioned  case,  it  is  obvious 


that  the  trajectory  obtained  from  multiple  shooting  was 
sub-optimal. 


Fig.  5.  Trajectories  under  ill-conditioned  case  with  selection 
of  IV  =  201,  IVa  =  30,  Ns  =  16.  The  trajectory  obtained 
from  local  pursuit  was  essentially  optimal,  while  the  one 
obtained  from  multiple  shooting  was  far  away  from  optimum. 


4-2.3  Large  number  of  segments 

When  N  >  600  the  orbit  transfer  problem  could  not  be 
solved  a  PC  with  1Gb  of  RAM  using  MS  3  .  On  the  other- 
hand,  PBMS’s  lower  memory  requirements  meant  that 
the  algorithm  was  able  to  operate  and  converge  to  the 
optimum.  Here  we  used  Na  =  60,  Ng  =  32;  the  final 
trajectory  is  shown  in  Fig.  6. 

5  Conclusions  and  ongoing  work 

This  paper  explored  the  use  of  a  bio-inspired  coopera¬ 
tive  strategy  termed  “local  pursuit”  for  solving  a  class 
of  numerical  optimal  control  problems.  We  discussed  a 
pursuit  algorithm  appropriate  for  problems  with  fixed 
final  time  and  partially-constrained  final  states.  The  al¬ 
gorithm  mimics  the  foraging  behavior  of  ant  colonies  and 
allows  a  collective  to  discover  optimal  controls,  starting 
from  an  initial  suboptimal  solution  and  using  only  local, 
pair-wise  interactions. 

We  proposed  combining  local  pursuit  with  multiple 
shooting  (MS)  as  a  way  to  overcome  some  of  the  com¬ 
putational  and  storage  limitations  of  the  latter  method. 


3  This  includes  the  memory  requirements  of  the  function 
used  to  compute  condition  numbers.  If  we  did  not  want  to 
record  condition  numbers  and  only  used  Gaussian  elimina¬ 
tion  method  to  solve  the  linear  system,  then  N  could  be 
increased  a  bit  further. 
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Fig.  6.  Spacecraft  trajectory  in  the  case  of  large  number  of 
segments,  N  =  601,  Na  =  60,  Ns  =  32. 


The  new  method,  termed  pursuit-based  multiple  shoot¬ 
ing  (PBMS)  involves  a  kind  of  “breaking  down”  of  the 
original  problem  down  to  smaller  segments,  to  be  op¬ 
timized  (using  MS)  by  a  group  of  (simulated)  agents. 
PBMS  has  slower  convergence  than  pure  MS  but  can  ex¬ 
ceed  the  limitations  of  MS  due  to  memory  requirements 
by  taking  advantage  of  the  cooperation  among  a  group, 
and  of  the  nonlinear  relationship  between  a  problem’s 
time  horizon  and  the  memory  required  to  solve  it  by 
MS.  Furthermore,  for  an  optimal  control  problem  with 
a  given  time  horizon,  PBMS  gives  one  the  option  to  in¬ 
crease  the  number  of  intermediate  points  along  the  solu¬ 
tion,  increasing  the  solution’s  accuracy  and  performing 
better  in  cases  where  MS  becomes  ill-conditioned.  These 
properties  of  PBMS,  and  a  comparison  with  MS  were 
illustrated  using  an  example  involving  orbit-transfer  of 
a  spacecraft. 


6  Acknowledgment 

The  authors  would  like  to  thank  Prof.  P.  S.  Krish- 
naprasad  and  Prof.  B.  Balachandran  for  their  helpful 
discussion  on  the  subject  of  numerical  computation. 


References 

[1]  J.T.  Betts.  Survey  of  numerical  methods  for  trajectory 
optimization.  Journal  of  Guidance,  Control  and  Dynamics , 
21  (2):  193-207,  Mar. -Apr.  1998. 

[2]  A.M.  Bruckstein.  Why  the  ant  trails  look  so  straight  and 
nice.  The  Mathematical  Intelligencer ,  15(2):59-62,  1993. 

[3]  A.M.  Bruckstein,  C.L.  Mallows,  and  I.  A.  Wagner. 
Probabilistic  pursuits  on  the  grid.  The  American 
Mathematical  Monthly ,  104(4) :323-343,  April  1997. 


[4]  S.  Camazine,  J.-L.  Deneubourg,  N.  R.  Franks,  J.  Sneyd, 
G.  Theraulaz,  and  E.  Bonabeau.  Self- Organization  in 
Biological  Systems.  Princeton  University  Press,  Princeton, 
New  Jersey  08540,  2001. 

[5]  M.  Dorigo,  V.  Maniezzo,  and  A.  Colorni.  Ant  systems: 
Optimization  by  a  colony  of  cooperating  agents.  IEEE 
Transactions  on  Systems,  Man  and  Cybernetics,  Part  B , 
26(1):29-41,  1996. 

[6]  D.M.  Gordon.  Ants  at  work.  The  Free  Press,  New  York,  1999. 

[7]  D.  Hristu-Varsakelis.  Robot  formations:  Learning  minimum- 
length  paths  on  uneven  terrain.  In  Proceedings  of  the  8th 
IEEE  Mediterranean  Conference  on  control  and  Automation, 
2000. 

[8]  D.  Hristu-Varsakelis  and  C.  Shao.  Biologically-inspired 
optimal  control:  learning  from  social  insects.  the 
International  Journal  of  Control,  77(18):  1545-1566,  Dec. 

2004. 

[9]  A.  Jadbabaie,  J.  Lin,  and  A.S.  Morse.  Coordination  of  groups 
of  mobile  autonomous  agents  using  nearest  neighbor  rules. 
IEEE  Transactions  on  Automatic  Control,  48(6),  June  2003. 

[10]  R.  Kress.  Numerical  Analysis.  Springer-Verlag  New  York 
Inc,  New  York,  1998. 

[11]  R.  Kurazume  and  S.  Hirose.  Study  on  cooperative  positioning 
system:  optimum  moving  strategies  for  cps-iii.  In  Proceeding 
of  the  1998  IEEE  International  Conference  in  Robotics  and 
Automation,  volume  4,  pages  2896-2903,  Leuven,  Belgium, 
May  1998. 

[12]  A.  Neumaier.  Solving  ill-conditioned  and  singular  linear 
systems:  a  tutorial  on  regularization.  SIAM  Review , 
40(3):636-666,  Sep.  1998. 

[13]  J.K.  Parrish  and  W.M.  Hammer.  Animal  groups  in  three 
dimensions.  Cambridge  University  Press,  Cambridge,  U.K, 
1997. 

[14]  K.  Passino,  M.  Polycarpou,  D.  Jacques,  M.  Pachter,  Y.  Liu, 
Y.  Yang,  M.  Flint,  and  M.  Baum.  Cooperative  control  for 
autonomous  air  vehicles.  In  R.  Murphey  and  P.M.  Pardalos, 
editors,  Cooperative  control  and  optimization,  pages  233-272. 
Kluwer  Academic  Publishers,  2002. 

[15]  S.I.  Roumeliotis  and  G.A.  Bekey.  Distributed  multi¬ 
robot  localization.  IEEE  Transactions  on  Robotics  and 
Automation,  18(5):781-795,  2002. 

[16]  C.  Shao  and  D.  Hristu-Varsakelis.  A  local  pursuit  strategy  for 
bio-inspired  optimal  control  with  partially-constrained  final 
state.  Technical  Report  TR  2005-76,  Institute  for  Systems 
Research,  University  of  Maryland,  College  Park,  MD  20742, 

2005. 

[17]  S.R.  Vadali  and  R.  Nah.  Fuel-optimal  planar  earth-mars 
trajectories  using  low-thrust  exhaust-modulated  propulsion. 
Journal  of  Guidance,  Control,  and  Dynamics,  23(3):476-482, 
May-Jun.  2000. 

[18]  I.A.  Wagner,  M.  Lindenbaum,  and  A.M.  Bruckstein. 
Distributed  covering  by  ant-robots  using  evaporating  traces. 
IEEE  Transactions  on  Robotics  and  Automation,  15(5) :918- 
933,  1999. 


10 


